home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / obsolete / pde_adi.pro < prev    next >
Text File  |  1997-07-08  |  7KB  |  250 lines
  1. ; $Id: pde_adi.pro,v 1.3 1997/01/15 04:02:19 ali Exp $
  2. ;
  3. ;  Copyright (c) 1992-1997, Research Systems Inc.  All rights
  4. ;  reserved. Unauthorized reproduction prohibited.
  5.  
  6. function BOUNDARY_LT,y,t
  7. ; left boundary condition as a function of y and t
  8.   bound=75.0
  9.   return,bound
  10. end
  11.  
  12. function BOUNDARY_RT,y,t
  13. ; right boundary condition as a function of y and t
  14.   bound=50.0
  15.   return,bound
  16. end
  17.  
  18. function BOUNDARY_BOT,x,t
  19. ; bottom boundary condition as a function of x and t
  20.   bound=0.0
  21.   return,bound
  22. end
  23.  
  24. function BOUNDARY_TOP,x,t
  25. ; top bondary condition as a function of x and t
  26.   bound=100.0
  27.   return,bound
  28. end
  29.  
  30. function NONHOMOGENEOUS,x,y,t
  31. ; the right hand side of the differential equation
  32. ; as a function of x, y, and t
  33.   nonhomo=0.0
  34.   return,nonhomo
  35. end
  36.  
  37. function INITIAL_COND,x,y
  38. ; the initial condition as a function of x and y
  39.   initial=0.0
  40.   return,initial
  41. end
  42.  
  43.  
  44. PRO PDE_ADI 
  45. ;+        
  46. ; NAME:
  47. ;    PDE_ADI
  48. ;
  49. ; PURPOSE:
  50. ;    This procedure computes the numerical solution of the 
  51. ;    initial-boundary value problem (parabolic partial  
  52. ;       differential equation) using the ADI method
  53. ;       (ALTERNATING DIRECTION IMPLICIT) in two space dimensions (x,y).
  54. ;      
  55. ; CATEGORY:
  56. ;    Numerical Analysis
  57. ;    
  58. ; CALLING SEQUENCE:
  59. ;    .run PDE_ADI     (compiles the procedure)
  60. ;    PDE_ADI        (runs the procedure)
  61. ;
  62. ; INPUTS:
  63. ;    There are no input parameters that can be specified from
  64. ;    the IDL command line. However, the following Internal
  65. ;    Parameters should be specified by editing the PDE_ADI.PRO
  66. ;    file before compiling the procedure:
  67. ;
  68. ;       DIFFUSIVITY ......................... DIFF
  69. ;       ENDPOINTS ........................... XMAX, YMAX
  70. ;       TIME STEP ........................... DELTA_TIME
  71. ;       NUMBER OF INTERIOR X-NODES .......... X_NODE_MAX
  72. ;       NUMBER OF INTERIOR Y-NODES .......... Y_NODE_MAX
  73. ;       NUMBER OF TIME STEPS ................ TIME_MAX
  74. ;
  75. ; OUTPUTS:
  76. ;    The numerical solution computed at each time step is stored
  77. ;    in the file PDE_ADI.DAT.
  78. ; 
  79. ; RESTRICTIONS:
  80. ;    All Internal Parameters must be positive scalars.
  81. ;
  82. ; PROCEDURE:
  83. ;    This procedure uses the ADI method to compute the  
  84. ;       temperature U at each node (M,N) in the rectangular  
  85. ;       grid defined by DELTA_X and DELTA_Y for each time 
  86. ;       step from time=0 to time=TIME_MAX with steps of 
  87. ;       DELTA_TIME.
  88. ;
  89. ; EXAMPLE:
  90. ;    This procedure solves the partial differential equation:
  91. ;       Ut-diff*(Uxx+Uyy)=NONHOMOGENEOUS(X,Y,T)
  92. ;       Subject to the initial and boundary conditions:
  93. ;       U(x,y,0)=INITIAL_COND(X,Y)
  94. ;       U(0,y,t)=BOUNDARY_LT(Y,T)   U(1,y,t)=BOUNDARY_RT(Y,T)
  95. ;       U(x,0,t)=BOUNDARY_BOT(X,T)  U(x,1,t)=BOUNDARY_TOP(X,T)
  96. ;       where:
  97. ;             Ut ....   first partial derivative of U with respect to t
  98. ;             Uxx ...  second partial derivative of U with respect to x
  99. ;          Uyy ...  second partial derivative of U with respect to y
  100. ;
  101. ;       The computational grid has X_NODE_MAX x Y_NODE_MAX interior mesh points
  102. ;       and X_NODE_MAX+2 x Y_NODE_MAX+2 total mesh points. 
  103. ;
  104. ; MODIFICATION HISTORY:
  105. ;    GGS; AUGUST 1992
  106. ;    Adapted from:
  107. ;        Graduate course work at the University of Colorado, Boulder
  108. ;    Published References:
  109. ;        Advanced Engineering Mathematics (sixth edition)
  110. ;        Erwin Kreyzig
  111. ;        Wiley & Sons Inc.
  112. ;
  113. ;-
  114.  
  115. ON_ERROR,2
  116.  
  117. ; start timer
  118. TIME_0=SYSTIME(1)
  119.  
  120. ; define internal parameters
  121. DIFF=1.0
  122. DELTA_TIME=.01
  123. XMAX=1.0
  124. YMAX=1.0
  125. X_NODE_MAX=8
  126. Y_NODE_MAX=8
  127. TIME_MAX=50
  128.  
  129. ; allocate arrays
  130. X=FLTARR(X_NODE_MAX+2)
  131. Y=FLTARR(Y_NODE_MAX+2)
  132. U=FLTARR(X_NODE_MAX+2,Y_NODE_MAX+2)
  133. V=FLTARR(X_NODE_MAX+2,Y_NODE_MAX+2)
  134. DIMENSION=max([X_NODE_MAX,Y_NODE_MAX])+1
  135. A=FLTARR(DIMENSION)
  136. B=FLTARR(DIMENSION)
  137. C=FLTARR(DIMENSION)
  138. D=FLTARR(DIMENSION)
  139.  
  140. ; define computational grid
  141. DELTA_X=XMAX/FLOAT(X_NODE_MAX+1)
  142. DELTA_Y=YMAX/FLOAT(Y_NODE_MAX+1)
  143.  
  144. ; define computational ratios
  145. X_RATIO=DIFF*DELTA_TIME/(DELTA_X^2.0)
  146. Y_RATIO=DIFF*DELTA_TIME/(DELTA_Y^2.0)
  147.  
  148. T=0.0
  149.  
  150. ; realize computational grid
  151. FOR M=1,X_NODE_MAX DO $
  152.   X(M)=X(M-1)+DELTA_X
  153. FOR N=1,Y_NODE_MAX DO $
  154.   Y(N)=Y(N-1)+DELTA_Y
  155. X(X_NODE_MAX+1)=XMAX
  156. Y(Y_NODE_MAX+1)=YMAX
  157.  
  158. ; initialize corners of computational grid 
  159. U(0,0)=0.5*(BOUNDARY_LT(Y(0),T)+BOUNDARY_BOT(X(0),T))
  160. U(0,Y_NODE_MAX+1)=0.5*(BOUNDARY_LT(Y(0),T)+BOUNDARY_TOP(X(0),T))
  161. U(X_NODE_MAX+1,Y_NODE_MAX+1)=0.5*(BOUNDARY_TOP(X(0),T)+BOUNDARY_RT(Y(0),T))
  162. U(X_NODE_MAX+1,0)=0.5*(BOUNDARY_RT(Y(0),T)+BOUNDARY_BOT(X(0),T))
  163.  
  164. ; compute initial values of U(M,N) on the interior of grid
  165. ; this is the temperature of the interior grid at TIME= zero
  166. FOR N=1,Y_NODE_MAX DO BEGIN
  167.   FOR M=1,X_NODE_MAX DO $
  168.     U(M,N)=INITIAL_COND(X(M),Y(N))
  169. ENDFOR
  170.  
  171. ;open data file: PDE_ADI.DAT
  172. OPENW,22,'PDE_ADI.DAT'
  173.  
  174. ; output solution at T=0.0
  175. PRINTF,22,FORMAT='(/"  Time = ",F5.2,/)',T
  176. PRINTF,22,'  NUMERICAL SOLUTION'
  177. FOR N=Y_NODE_MAX+1,0,-1 DO $
  178.   PRINTF,22,FORMAT='(20(F7.2,1x))',U(0:X_NODE_MAX+1,N)
  179.  
  180. ; start time steps
  181. FOR J=1,TIME_MAX DO BEGIN
  182.  
  183. ; start y-direction computation
  184.   FOR N=1,Y_NODE_MAX DO BEGIN
  185.     FOR M=1,X_NODE_MAX DO BEGIN
  186.       A(M)=-X_RATIO/2.0
  187.       B(M)=1.0+X_RATIO
  188.       C(M)=-X_RATIO/2.0
  189.       D(M)=0.5*Y_RATIO*U(M,N-1)+(1.0-Y_RATIO)*U(M,N)+0.5*Y_RATIO*U(M,N+1)
  190.       D(M)=D(M)+0.5*DELTA_TIME*(NONHOMOGENEOUS(X(M),Y(N),T) $
  191.                +NONHOMOGENEOUS(X(M),Y(N),T+DELTA_TIME))
  192.     ENDFOR  
  193.     U(0,N)=BOUNDARY_LT(Y(N),T+DELTA_TIME)
  194.     V(0,N)=0.25*Y_RATIO*(BOUNDARY_LT(Y(N-1),T)+BOUNDARY_LT(Y(N+1),T)) $
  195.                +(0.5-0.25*Y_RATIO)*BOUNDARY_LT(Y(N),T) $
  196.                -0.25*Y_RATIO*(BOUNDARY_LT(Y(N-1),T+DELTA_TIME) $
  197.                +BOUNDARY_LT(Y(N+1),T+DELTA_TIME)) $
  198.                +(0.5+0.25*Y_RATIO)*BOUNDARY_LT(Y(N),T+DELTA_TIME)
  199.     D(1)=D(1)+0.5*X_RATIO*V(0,N)
  200.     U(X_NODE_MAX+1,N)=BOUNDARY_RT(Y(N),T+DELTA_TIME)
  201.     V(X_NODE_MAX+1,N)=0.25*Y_RATIO*(BOUNDARY_RT(Y(N-1),T)+BOUNDARY_RT(Y(N+1),T)) $
  202.                           +(0.5-0.25*Y_RATIO)*BOUNDARY_RT(Y(N),T) $
  203.                           -0.25*Y_RATIO*(BOUNDARY_RT(Y(N-1),T+DELTA_TIME) $
  204.                           +BOUNDARY_RT(Y(N+1),T+DELTA_TIME)) $
  205.                           +(0.5+0.25*Y_RATIO)*BOUNDARY_RT(Y(N),T+DELTA_TIME)
  206.     D(X_NODE_MAX)=D(X_NODE_MAX)+0.5*X_RATIO*V(X_NODE_MAX+1,N)
  207.     TRIDAG,A(1:X_NODE_MAX),B(1:X_NODE_MAX),C(1:X_NODE_MAX),D(1:X_NODE_MAX),WEIGHT
  208.     FOR M=1,X_NODE_MAX DO $
  209.       V(M,N)=WEIGHT(M-1)
  210.   ENDFOR
  211.  
  212. ; start x-direction computation
  213.   FOR M=1,X_NODE_MAX DO BEGIN
  214.     FOR N=1,Y_NODE_MAX DO BEGIN
  215.       A(N)=-Y_RATIO/2.0
  216.       B(N)=1.0+Y_RATIO
  217.       C(N)=-Y_RATIO/2.0
  218.       D(N)=0.5*X_RATIO*V(M-1,N)+(1.0-X_RATIO)*V(M,N)+0.5*X_RATIO*V(M+1,N)
  219.       D(N)=D(N)+0.5*DELTA_TIME*(NONHOMOGENEOUS(X(M),Y(N),T) $
  220.                +NONHOMOGENEOUS(X(M),Y(N),T+DELTA_TIME))
  221.     ENDFOR
  222.     U(M,0)=BOUNDARY_BOT(X(M),T+DELTA_TIME)
  223.     D(1)=D(1)+0.5*Y_RATIO*U(M,0)
  224.     U(M,Y_NODE_MAX+1)=BOUNDARY_TOP(X(M),T+DELTA_TIME)
  225.     D(Y_NODE_MAX)=D(Y_NODE_MAX)+0.5*Y_RATIO*U(M,Y_NODE_MAX+1)
  226.     TRIDAG,A(1:Y_NODE_MAX),B(1:Y_NODE_MAX),C(1:Y_NODE_MAX),D(1:Y_NODE_MAX),WEIGHT
  227.     FOR N=1,Y_NODE_MAX DO $
  228.       U(M,N)=WEIGHT(N-1)
  229.   ENDFOR
  230.  
  231. ; let user know how the algorithm is progressing
  232.    PRINT,' ... time step complete     REAL TIME = ',T, ' seconds'
  233.  
  234. ; advance to next time step
  235.     T=T+DELTA_TIME
  236.  
  237. ; output solution at time step T
  238.     PRINTF,22,FORMAT='(/"  Time = ",F5.2,/)',T
  239.     PRINTF,22,'  NUMERICAL SOLUTION'
  240.     FOR N=Y_NODE_MAX+1,0,-1 DO $
  241.       PRINTF,22,FORMAT='(20(F7.2,1x))',U(0:X_NODE_MAX+1,N)
  242. ENDFOR    
  243.  
  244. ; close data file: PDE_ADI.DAT
  245. CLOSE,22
  246.  
  247. PRINT,'             program execution in ...',SYSTIME(1)-TIME_0,' seconds'
  248.  
  249. END
  250.